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COMPUTER PROGRAM FOR FLAT SECTOR 
THRUST BEARING PERFORMANCE 
by Alden F. Presler and Izhak Etsion 
Lewis Rv'isearch Center 

SUMMARY 

A versatile computer program is presented which achieves a rapid, numer- 
ical solution of the Reynolds equation for a flat sector thrust pad bearing with 
either compressible or liquid lubricants. Program input includes a range in 
values of the geometric and operating parameters of the sector bearing. Per- 
formance characteristics are obtained from the calculated bearing pressure 
distribution. These are the load capacity, center-of-pressure coordinates, fric- 
tional energy dissipation, and flow rates of liquid lubricant across the bearing 
edges. Two sample problems are described and solved, one each for gas and 
liquid lubricants. 


INTRODUCTION 

The computer program presented in this report calculates the performance 
characteristics of both gas and liquid lubricated flat sector pad thrust bearings. 

A typical pad configuration is shown in figure 1. The bearing consists of 
several pads each of which has an extent angle (3 and inner and outer radii r^ 
and r^, respectively. Each pad assumes both pitch and roll with respect to the 
rotating runner to provide a (generally) converging film thickness in the direc- 
tion of rotation. 

In a recent paper (ref. 1) it was shown that any pitch and roll of a sector 
shaped pad about a certain point can be transformed to a corresponding pure 
pitch about a certain radial line. This can be imderstood from figure 1 by visual- 
izing a plane parallel to the runner that goes through the origin of the sector 
(point 0 in fig. 1). The radial pivot line is the intersection between this parallel 
plane and the plane of the tilted sector, and it can be either within or outside the 
sector boundaries. 

Based on this observation the flat sector thrust pad has been analyzed for 
both compressible and incompressible lubricants in references 1 and 2. 

The objects of this report are to: (1) describe the numerical analysis and 
solution of the basic Reynolds' equation, (2) document the resulting computer 
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programs FSTBPl and FSTBP2 for compressible (gas) and incompressible 
(liquid) lubricants, respectively, and (3) serve as a users guide for these two 
programs. 

Input data for the programs describe the physical characteristics and the 
geometry of the sector pad. These are inner- to- outer radius ratio, pad angle 
extent, pad pivot angle, the ratio of pad pitch to minimum film thickness, and 
bearing compressibility number (this is not specified for liquid films). 

The computed results include the pad thrust loading, the frictional power 
loss coefficient, and the center- of-pressure coordinates. Additional calculated 
results are flow leakage from the downstream edges for the liquid film case. 

(For a complete bearing the results of the several pads are added. ) 

This report also contains six appendices which give complete detsils of the 
numerical methods of solution, two FORTRAN listings of the computer programs, 
and two sample problems with output listings. 


STATEMENT OF THE PROBLEM 

It is required to develop a computer program that numerically solves the 
lubrication boundary value problem defined physically by the Reynolds equation 
over a tilted sector pad. 

The Reynolds equations are essentially diffusion equations for the lubricant 
film pressure. In cylindrical coordinates the equation for the isothermal com- 
pressible case with density proportional to pressure is 

S, ^] + l ±[ 2 ^ = 6ru> (1) 

8r \ ju dr J r 90 \ fj, ddj 90 

and for the incompressible lubricant film 

( 2 ) 

9r \ p drj r 90\/it 90y 90 

The lubricant film thickness is expressed in terms of the independent variables 
r and 0. By considering the clearance h^ along the pivot line as a reference, 
the film thickness at any point (r, 0) is given by 

h = ho + yr sin(0p - 0) (3) 

where y is the amount of tilt or pitch about this line. All of the symbols used 
in these and the following equations are defined in appendix A. 
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The normalized (dimensionless) form of these equations are 



for compressible films, and 

~ ^ +1 

3R \ aR/ R ae \ 

for an incompressible lubricant. From equation (3), the film thickness becomes 

H = 1 + eR sin(6»p - 0) 

Normalization of the variables leading to equations (4) to (6) is described in 
references 1 and 2. 

The pressure distribution over the pad area is obtained by numerically solv- 
ing equation (4) or (5) under boundary conditions which are defined in the next 
section. Further sector pad calculations are based on this pressure distribution. 
The following bearing performance characteristics are then calculated; 

(1) Pad load W 

(2) Center- of- pressure coordinates R^p, 9^^ 

(3) Power loss coefficient (normalized coefficient of friction) F/ W 

(4) Volumetric lubricant flow rates across the sector pad edges (liquids only) 

•lie- V- «so- '•si 

Inputs are in vector arrays with a range of design parameter values: 

(1) Pad dimensions - 

(a) Inner radius ratio R. = r./r^ >0 

(b) Sector angle p in degrees 

(2) Pivot line angle ratio 0p//3 

(3) Compressibility factor (also called ’’bearing number”) A, for the gas 

film case 

(4) Ratio of pad slope to minimum pad- runner clearance e/Hg > 0 

METHOD OF SOLUTION 

The two forms of the Reynolds equation are first transformed to the following 
boundary value problems, which are then solved numerically for specific boundary 
conditions. 


30 J 89 
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Find a function u(R, 0) that satisfies the equation 

L(u) = f(R, 0) 

on the domain 

= (0, R 0 < 0 ^ j3; R. s R ^ 1; /3 - - < 0 <-] 

\ ^ 2 P 2/ 

In the compressible ease, L is the nonlinear operator 






3R^ R"' 90" 


2 ^.2 \9R R/ 9R \ 00 R 


ah; 






H 




00 


and 


u = Q = (PH)"= 
f = 0 

u = on boundary 8^, (P - 1) 


( 7 ) 


(8) 


/ (9) 




In the incompressible case, L is the linear operator 


and 


L=-ii+i-^ + 

0R^ R^ 00^ 


In RH^\ In RH^|-^ 

,0R / 0R ^j^2 00 y00 




u = P 



u = 0 on boundary 0/&, (P = 0) 


(10) 


Note again that on the pad boundary P = 1 (compressible) and P = 0 (incom- 
pressible). On the boundary 8a'^ the pressure is fixed at the ambient value, 

P = Po> ^*1 liquid lubricant cases. But the definition of dimension- 

A ' • ■ 

less pressure differs-. 


P = - 2 - (compressible) 
P, 


P = 


a 

P' P. 




( 11 ) 


K 


( in compre s sible) 


where 


K = 


6jjLU)r 


Numerical solution of equation (9) or (10) is described in appendix B. These 
boundary value equations are approximated by finite difference equations on a 
polar mesh over the sector pad area as indicated in figure 2. In order to adopt 
Simpson Rule integration formula, the R and 6 intervals are divided into even 
numbers of increments. 

As indicated in the previous section, the solution of the equations is in the 
form of the film pressure distribution over the pad area. From this distribution 
the thrust pad performance characteristics are obtained by Simpson Rule integra- 
tion of the following e3q>ressions: 

Normalized load capacity 


^ ^ = f (P - 1) R d0 dR (compressible) 


W = 


„ 2 'R. 0 

Pr. 1 

*^a o 


(12a) 


or 


W = — ^ = f f ^ PR d0 dR (incompressible) 
<> J-R ■'n 


Kr2 "«i "0 


(12b) 


Unit load 


W 

UL = (compressible) 

V 


W 

UL = — — (incompressible) 
KA 


2W 


/3(l - R?) 


(13) 


Where the W is calculated from equation (12a) or (12b) depending on fluid type. 
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Center- of-pressure radial coordinate 


^cp 


. 


, ( P r i ) R: .. dR . d^ (compressible) 
W 


(14a) 


or 


R 


2 

PR^ dR d0 


Cp 


W 


(incompr es sible) 


(14b) 


>0 UR. 


Center- of- pres sure angle coordinate 


= sin” ^ 
cp 


/? pi 


^i 


(P - pR^^ dR sin 6> dg 

R„ W 
cp 


(compressible) (1 5a) 


or 


cp 


^ ^ ‘ Pr 2 dR sta e de \ ,i„eompresslble) 


•0 \/R. 


cp 


(15b) 


Center- of-pressure distance from pivot line 


= R^„ sin(0„^ - 0 \ 
cp cp . ' cp p' 


(16) 


Power loss coefficient 


H =_1 

W 6W 


'/? ri 


'0 ^R. 


+ 3R-ii- ^\dR de 


ABT , H 


H/H2 Hg de j 


(17) 


Equation (17) holds for both compressible fluids (A variable) and incom- 
pressible fluids. In the incompressible case, the program puts A = 1. 
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Volumetric flow rates across the pad edges are: 
Trailing edge; 



Leading edge: 





Flow is defined as positive in the direction of increasing R and 6 . 

FORTRAN PROGRAM 
General Description 

The foregoing analysis has resulted in two thrust pad computer programs; 
FSTBPl for compressible fl'oid films, and FSTBP2 for incompressible lubricant 
films. The FORTRAN listing for each program is given in appendix C, and flow 
chart diagrams are presented in appendix D. The dictionary of the FORTRAN 
symbols used in the computer programs is appendix E. 


-/ 


r 




'1 


I 


T 


J 
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Both programs have identical structure in the number and function of the 
subprograms, and in the general format of the input data. 

The first flow diagram in appendix D, that of the supervisory module MAIN2, 
presents a compact overview of the logical sequence which the computer programs 
follow in producing pad performance characteristics from input operating condi- 
tions. 

The first part of the computer program, which includes the first four sub- 
routines called by the module MAIN2, accomplishes the numerical solution of the 
Reynolds equation using a Gauss-Seidel iterative method with under- or over- 
relaxation. Upon convergence of the iterations, the calculated film pressure dis- 
tribution under the sector pad is passed to the second part of the program which 
performs the numerical double integrations on equations (12) to (17) resulting in 
the bearing pad performance characteristics. 

A brief explanation of the function of each of the program modules follows: 

(1) MAIN2 is the executive routine for processing multiple cases, and has 
primary control of logical flow throughout the complete program. The executive 
routine also controls the printing of the res^ilts. 

(2) EUCLID creates the mesh AR, A0,' and converts all angles to radians. 

(3) XBEGN2 tests the given angular position of the pad pivot line against the 
pad coordinates, and then calculates the minimum film thickness between pad and 
runner, the pad slope, and the maximum- to- minimum film thickness ratio. 

(4) COEFF generates the values of the nodal coefficients for the finite differ- 
ence representation of the Reynolds equation. It also provides the initial values 
of the dependent variable for the first iteration of the Gauss-Seidel process. 

(5) RELAX is the basic working routine for solving the Reynolds equation by 
Gauss-Seidel iterative method with a choice of relaxation parameter, 

(6) TABULT is an integration subroutine completing the calculation of the 
pad loading, center-of-pressure coordinates, friction power loss, and lubricant 
flow rates. 

(7) RSIMP uses the Simpson Rule to Integrate tabulated functions along the 
radii at each angle mesh position. 

No special effort was made to determine the optimum relaxation factor for 
the program with given input data. Initially the sequence for = k/4, where 
k was incremented from 1 to 7, was tested for iteration efficiency. It was soon 
determined that = 1. 25 was generally superior to other values tried, and 
that value was used for most calculations. As compressibility number increased 
to 100 the optimum value of 9 . decreased to 0. 75. 


/ 
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Using the Program 

The program always starts by reading in a three- card input deck. The con- 
tents of these cards is now described (all symbols are repeated in appendix E). 

Card 1 (FORMAT 316, F8. 2, 2E8, 1) 

NR number of radial increments (even) 

NA number of angular increments (even) 

INTRMX maximum number of iterations allowed if convergence is not 
achieved in the solution. 

OMEGA relaxation parameter in the solution procedure. 

HALT cutoff value for minimum film thickness. There is an error return 

and a warning printout if the calculated HMIN is less than HALT. 

RESIDL convergence criterion; test on maximum change in dependent vari- 
able between iterations in subroutine RELAX. 

Card 2 (FORMAT 4L6) 

Values on this card control logical switches. 

DEBUG when TRUE the values of the dependent variable are printed for all 
mesh points for each iteration. To avoid printout use FALSE. 

TABOUT when TRUE the film pressure array over all mesh points is printed 
after convergence of the iteration routine. The value FALSE 
stops this printout. 

OLDQ a TRUE value saves the dependent variable array after convergence 

of the iterative procedure so that this array can be used as the 
starting value for a new iterative calculation when multiple input 
cases are being calculated. A value of FALSE defines initial 
pressure as ambient over pad area. 

VARGRD a TRUE value creates a variable mesh grid over the sector pad 
(fig. 2). Only needed for liquid lubricant calculations. 

Card 3 (FORMAT 5110) 

This card sets the size indices for parameter arrays which are read in imme- 
diately after this card through namelist '’VARBLE. *» 
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NTH-ETA 

NRATIO 

NUMRI 

NUMBET 

NUMLMB 


quantity of 6 ^/^ values 
number of e/Hg values 
selection of RI ratios 
number of (3 values 

how many A; can be left blank for liquid calculations. 


The parameter data whose array sizes are defined by the values on card 3 
are read into the MAIN2 subprogram throu^ the namelist " VARBLE'* input, 
utilizing as many cards as is necessary. These are without format description, 
and are now listed. 

TRATIO (I) 

1=1, NTHETA 

ERATIO (I) 

1=1, NRATIO 

VRI (J) 

J=l, NUMRI 

VBETA (J) 

J=l, NUMBET 

VLMBDA (J) 

J=l, NUMLMB 

Documented description of the data input decks are provided at the end of 
each of the two programs FSTBPl and FSTBP2 in appendix C. 

CONCLUDING REMARKS 

We have described a numerical method and computer program for solving 
two forms of the Reynolds equation within a circular sector region. These two 
forms of the equation are for compressible and incompressible fluid films. 

The numerical method uses a curvilinear cell at each mesh point to derive 
the f inite difference analog to the Reynolds equation. This represents a system 
of nonlinear equations with prescribed constant boundary values. Two computer 
programs were developed to solve the finite difference systems representing the 
compressible and incompressible fluid cases. The programs use a Gauss-Seidel 
iterative method with relaxation capability. 


array of pivot line angles 0 //3 

P 

array of input parameters yro/h2 = 
array of inner to outer radius ratios, r./ r^ 
array for pad angle size /3 


for compressible films, the array of bearing numbers A 
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liiput flexibility was built into the programs to allow a wide variation in sec- 
tor geometry and operating conditions. This allows the program user to quickly 
evaluate alternative thrust bearing design configurations in order to find the 
optimum cases. 

The programs are easy to use, and have enough generality to be used as the 
"black box" in situations where steady state solutions to the Reynolds equation 
are required for circular sector regions. 

A rather coarse mesh of about 100 points is sufficient for the bearing load 
and moment calculations (8 radial divisions and 12 angular divisions). Edge leak- 
age calculations which utilize pressure gradients at the boundaries require a finer 
mesh (14 radial and 20 angular divisions). 

Each case required at most 40 to 50 iterations to satisfy the convergence 
criterion, and averaged 1 to 2 seconds total time on the UNIVAC 1100/42 
computer. 
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A 


8h 

F 

F 


F/W 

f 

H 


Hi 

«2 

h 

hi, 

K 


L 


P 

P 


Pa 

Q 

q 

R 


APPENDIX A 

MATHEMATICAL SYMBOLS 

pad area, /3(rQ-|r?y2 

closed domain oyer pad area 
boundary of pad domain 
friction power loss 

nondimensional px)wer loss, F/p^cohgr^ (incompressible) 

power loss coefficient, F/Wcohg 
forcing function in eq. (7) 
nondimensional film thickness, h./h 

maximum film thickness ratio, h^^/h^ 

minimum film thickness ratio, hg/h^ 

film thickness 

constant film thickness along the pivot line 
maximum and minimum film thickness, respectively 
6f.cr2/h2 

operator in Reynolds equation (eq. (7)) 

nondimensional pressure, p/p (compressible); or (p - p„)/K 

A O. 

(incompressible) 
pressure 
ambient pressure 
(PH)2 

volumetric edge flow rates, eqs. (18) to (21) 
nondimensional radius, r/r^ 
nondimensional radius to center of pressure 
nondimensional inner radius, r./r^ 



radial coordinate 
pad inner radius 

pad outer radius 

dependent variable at node (i, j) in difference eq, (E-3) 

unit load (eq. (13) = 2w/^^l - R? 

dependent variable in Rejuiolds equation (eq. (7)) 
pad load capacity 

nondinaensional load, w/p^r^ (compressible); or w/xr^ (incom- 
pressible) 

tilt about a tangent line (roll) 
tilt about a radial line (pitch) 
angular extent of pad 
pitch about pivot line 
finite difference increment 
tilt parameter, 

clearance parameter, 'Y'^q/^2 

angular coordinate 

2 / 2 

bearing compressibility number, G^wr^/p^hg 

relaxation factor 
bearing shaft speed 
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APPENDIX B 


NUMERICAL ANALYSIS 


Numerical solutions of the Reynolds equations (7) to (10) are obtained by 
discretizing the partial differential equations and solving the resulting set of 
difference equations by an iterative technique over the closed region -V of the 
sector pad area. 

The sector pad is partitioned into a polar mesh (fig. 2(a)) with an even num- 
ber of increments in both coordinates, resulting in an odd number of nodes. A 
variable mesh capability is provided for the liquid lubricant cases. The finest 
mesh starts at the pad boundary (fig. 2(b)). The enlarged view of the central 
difference mesh is shown in figure 3. 

In central difference formulations, the dependent variable U. . atnode (i, j) 
is a function of U at the four surrounding nodes. The central difference oper- 
ators for the Reynolds equations (7) to (10) are 




where 


aR'^ 


(AR)^ 




d6^ ~ 

nfiAef 

dV 

3 " ^i-1. 3 

8R 

2 AR 



39 

2A9 


A9=-^ I 


NA 1 


AR = 


1 - RI 
NR 


(B-1) 


(B-2) 


The errors in the finite difference operators (B-1) are all of the order (AR)" 
and (A^)^. 
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When these finite differences are substituted for the partial derivatives in 
the L- operator of equation (7), and the coefficient functions of the first order 
derivatives are expanded, the solution for U.j in terms of the LT s at the four 
surrounding nodes is obtained; 


, .+C..-U. +D..*U. . , +F.. 

u . = -U u IzL-J 1 ] 1.3+1 13 1,3-1 1 ] 

E.. 

ij 


(B-3) 


where the coefficients are given in table B-1. 

The solution of equation (B-3) proceeds from an initial estimate of the ujP 
over the mesh 1 :< i :< NR +1, 1 s j < NA + 1. It is always safe, if not efficient, 
to start with ambient pressure across the film mesh, that is. 




= Hjj (compressible) \ 
U-P = 0 (incompressible) 

J 


(B-4) 


Subsequent calculations of the U^j continue from the preceding values of the 
estimates using equation (B-3). In the GausS- Seidel iterative proce- 
dure, the are immediately replaced in the storage array by upon 

calculation, which modifies equation (B-3) to 


(n+1) ^ ^i+L 
ij 


+ B,,U) 


(n+1) 


+ c.,uy 


(n) 


J 12 ^-l.J 13 hJ+1 


+ D.. 
LL 


u(n+l) 


+ F.. 
11 


E.. 

13 


(B-5) 


over the domain 2 i < NR, 2 j NA. 

It is noted from table B-1 that for the compressible films the coefficients 

C - . and D-. contain the nonlinear term /u.. by way of G... This creates no 
J J V J 

problem as the preceding iterate is always used in the G^ for equation 

(B-5). 

Relaxation of the basic iterative numerical procedure is used in this report 
to hasten convergence of the calculations. The relaxation method modifies 
equation (B- 5) to the form 
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u>. 

1] 


(n+1) ^ 


E.. 

13 


(B^6) 


where the relaxation factor Q iS between 0 and 2. For « = 1, the regular 

Gauss- Seidel form of equation (B-5) results. 

Between successive iterations the change in the solution of equation (B-6) is 

denoted by 


TEST = 


:jj(n+l) _ 

ij 13 


(B^7) 


The largest value of TEST in one complete solution of equation (B-6) across 
the domain is stored as STRERR. At the completion of each solution on 
the maximum deviation or error between iterations is tested against a conver- 
gence criterion called RESIDE. When STRERR < RESIDL the numerical solution 

of equation (B-6) is defined as accomplished. 

The pressure distribution over the pad domain is determined from the 
last solution U.j after convergence using the definitions in equations (9) and (10); 








(compressible) 


(B-8) 


p„ = U- (incompressible) ^ 


The bearing performance characteristics W, R^p, and 0^^ are expressed 

by equations (11), (13), and (14), respectively. The numerical integrations are 
accomplishedusingSimpSotf s 1/3 Rule, and are in two steps. The first carries 
out the radial integration 


XF, 




K 


3 

NR/2 

i=l 


4(P2j j - l)R^i + 2(P2i+i_ j - DRfi+1 


K = 1, 2 


(B^9) 



The circumferential integration step completes the calculation: 


NA/2 

^1, K ■ ^NA+1, K K ^^2i+l, 

i=l 

where 

W=INTk=1 

and 



«cp = ^%=2 


(B-10) 


The calculation requires first 


AA 


A0 


NA/2 


XF. 


NA+1, 2 ^ 2 ®2i ^^^2i+l, 2 ®2i+l 

i=l 


(B^ll) 


and finally 


= sin 
cp 


AA 

(W.R^p) 


(&-12) 


The terms (P.. - 1) in equation (B-9) are correct only for the compressible 
case, and represent gage pressures. These terms must be just P.j for the in- 
compressible solution as is derived from the second expression in equation (B-8). 

The power loss coefficient, equation (17), and the leakage integrals, equa- 
tions (18) to (21), require the evaluation of the first derivatives of pressure along 
the domain boundary In order to achieve the same order of accuracy for the 

finite difference approximations on the boundaries as for the central difference 
expressions in equations (B-l), the following forward difference operation at the 
leading edge and inside radius are used; 


(^\ ~ 1 2 ~ ^i. 3 

Vae/g^Q 2A0 


(B-13) 
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F) 

V8R/ 


'e=e. 

R=R. 


P3.J 

2 AR 


(B-14) 


At the trailing edge, and along the outside radius of the pad, the backward differ- 
ence operators are used; 


^ ^i. NA-1 ~ ^^i. NA ^^i. NA+1 


dd I 


e=jS 

R.<R<1 


2 Ae 


(B-15) 


'^\ ^ %R-1. 1 ~ '^^R. 1 ^%R+1. 1 


8R/ 


R^l 


2 AR 


(&-16) 


As with the central difference operators (B-1), the forward and backward differ- 
ence operators (B-13) to (B-16) are of the order (AR)^ and (Afl)^ in error. 

Since the pressures in the gradient expressions can be gage pressures, and 
since gage = 0 on 3^, equations (B-13) to (B-16) can be further simplified 
to 


8P \ ^ ^^i. 2 ~ ^i. 3 

2A6 

R.<R<1 


8P\ ^ ^^2,j- P 3 J 


^^le=e. 

R=«i 


2 AR 


307, 


e=(3 

R.<R<1 


NA-1 ~ ^^i. NA 
2 A0 


- %R-l. i ~ ^%R.1 


3RA 


0 = 0 . 

R=1 


2 AR 
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TABLE B-1. - COEFFICIENTS FOR CENTRAL DIFFERENCE 
APPROXIMATION OF REYNOLDS EQUATION (EQ. (B-3)) 

Incompressible films 


Coeffi- 

cient 










Compressible films 

Uij = 


u. . = p. . 

13 1] 


“f • 


4- 


(AR)2 2.R.H...AR 


j (AR)^ ^i’ 




(AR)^ 2-R.H.^.-AR 


4- 


(AR)^ 2-R.-AR 




+ 


Rj(A0)^ 


I 


j r 2(A0)2 2.R.H...Ae 


e-cos(0^ - 0.) I 1 

2.R.Hy.A0 i 


e-cos(0 -0.) i 


; Rf(A0)^ 


3- e- cos(0p - 0^) 
2-R.H.j-A0 


3*e*cos(0p- 0^) 
2-R.H.yA0 


E 


ij 


1 1 . 

9* 

1 1 

_(AR)^ R?(A0)^_ 


_(AR)^ R?(A0)^_ 


^ij 


«ij 



Ry h|- e- cos(a - e.) 


H® 

"ij 




1 + €• Rj- sin(0p - 0^) 


1 + e- Rj- sin(0p- 6.) 
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APPENDIX C 


FORTRAN PROGRAMS 
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FSTBPl 

Flat Sector Pad Thrust Bearing Program Number 1 gas lubricant. 


>TPFf (0) .MAIN2 
C 

C**MAIN EXECUTIVE FOR GAUSS-SEIDEL ITERATION 
C 

REAL LANBDA, ONEGA 

LOGICAL DEBUG, TABOUT,OLDQ,VARGRD 

C 

COMNON/BLOGIC/DFBUG , TABOU T , OLOC , V ARGRD 
COMMON/GECM/ANGR , THPR,DR,DA 
DIMENSION EXIT(S) ,TRATIO( 15 ) ,ERATI0(26) 
DIMENSION VRI(10),VSETA(10),VLMBDA(1D) 
NAMELIST/VARBLE/TRATIOjERATIO, VRI ,VBETA, VLMBDA 


C 

c 


e 

10 




FORMAK 3I6,F8.2,2E8.1 ) 

FORMAK^Le) 

FORMAK 5I1U) 

FORMAK 2X , 12HN0 . OF R0WS=,I3,16H NO. OF C0LSr,I3,26H MAX NO. 
1 OF ITERATI0NS=,I14,25H RELAXATION PARAME TERr , FH . 2 , 19H RESIDU 
2AL ERROR:, F10.7/2X,36HSMALLEST ALLOWED FILM THICKNESS HALT r, 
3610.4//) 

FORMATi 2X,39HTHE MINIMUM FILM THICKNESS IS LESS TH AN ,G 1 0 . 4 / / ) 


READ (5, 1) NR,NA , ITERMX,0MEGA , HALT, RESIDE 
READ (5, 2) DEBUG, TABOU T,CLDQ, V ARGRD 
READ (5, 3)NTHETA ,NRATIO,NUMRI,NUMB£T,NUMLMB 
READ (5, VARBLE ) 


WRITE(6 ,4 )NR,NA,ITERMX,OMEGA,RESIDL,HALT 


WRITE(6 ,900) 

IPRINT=D 

DO 500 NRI=1,NUMRI 
RI=VRI(NRI) 

DO 4C0 NBETA=1 ,NUMBET 
BETA=VBETA(NBETA ) 

DO 300 NLMBD=1,NUMLMB 
LAMBDA:VLMBDA(NLMBO) 


S HEADING PRINTOUT. . 

9 INITIALIZE PRINT SWITCH... 

9 INCREMENT INNER RADIUS VALUES.. 
9 INCREMENT BETA VALUES. .DEGREES. 
9 INCREMENT BEARING NUMBER VALUES 


DO 200 NANGLE:l,NTHETA 
XXTHTP:TRATI0(NANGLE ) 

CALL EUCLID(NR,NA,PI,BETA,XXTHTP) 
KOUNT=D 


9 INCREMENT THETAP/BETA RATIO 

9 RETURNS GEOMETRIC PARAMS. 


IF(IPRINT.GT.D)G0 TO. 8 
WRITE (6, 910) 

GO TO 10 ' 

WRITE(6,915) 9PRINT START TOP OF PAGE.. 

WRITE(6,920) 

WRITE (6,930 >RI ,BET A, LAMBDA, XXTHTP 
WRITE! 6,94D) 

WRITE(6,950) 


DO 100 NFILM:1,NRATI0 


9 INCREMENT EPS/H2 RATIO 


o o o 


ERAT=ERATIO(NFILM> 
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C 

C 


C 

12 

C 

15 


C 


i: 

21 


C 

22 


C 


CALL BEGIN2(EPS,HALTfHMINtHRAT|ERAT,RI »I12> SFILM THICKNESS CALCS. 
60 TO 15 

WRITE(6,5)HMIN 9 ERROR EXIT.. 

60 TO IDO 

HLMBDA=L AMBDA«HMIN**2 
IF<K0UNT.6T.0)G0 TO 21 

CALL ARRAYS(NR»NA,EPS,KOUNT,RI,HLMBDA ) SFINITEDIFF. COEFFS. 

CALL RELAXfNR, NA, OMEGA, ITERMX, RESIDL ) SSOLVE REYNOLDS EQUATION, 

CALL TABULT(NR,NA,RI , EXIT, LAMBDA, HMIN ) SINTE6RALS OF PRESSURE 

aOVER PAD AREA., 


GO TO 22 

CALL RARRAY»EPS,KOUNT,HLMBDA ) SENTRY TO SUBROUTINE ARRAYS. 

CALL RRELAX SDITTO FOR THE 6AUSS-SEIDEL ROUTINE. 

CALL RTAB<EXIT, LAMBDA, HMIN» aENTRY TO INTEGRAL CALCS,. 

RR=EXIT(1) 

AArEXITI2) 

WW=EXITi3) 

FF=EXITC4) 

WUNIT=EXIT(5) 

XCP=RR*SIN(ANGB*AA-THPR1 

HCP=1.-EPS*XCP 

H2GP-HMIN/HCP 

HINVRSri ,/HMlN 

R3=ERAT 

R7=FF/WW 


WRITE(6,960)HRAT,HINVRS,R3,WW,UUNIT,FF,R7,H2CP,RR,AA,XCP 
C 

K0UNT=K0UNT^1 
IDO CONTINUE 

IPRINT=IPRINT+1 
200 CONTINUE 
C 

300 CONTINUE • 

.400 CONTINUE , ' 

“500 CONTINUE 
C 

900 FORMAK 1H1//40X,40HNUMERICAL RESULTS - GAS BEARI NG A N ALYS I S / 40X , 
128HLEWIS RESEARCH CENTER CNASA)////) 

910 FORMAT! SOX, 17HSYSTEM PAR AME TERS // 1 5X , 4H ( P 1 ) , 25X , 4 H ( P2 ) , 28 X , 4 H ( P3 ) , 
127X, 4H( P4 )//) 

915 FORMAT! 1H1//50X,17HSYSTEM P AR AMETERS//1 5X ,4 H !P1 » , 25X , 4H ! P2 ) , 28X , 
14H!P3) ,27X,4H!P4I// ) 

920 FORMAT! 10X,17HINNER/0UTER RADII ,14X,17HBEARING PAD ANGLE , 14X , 14HBE 
1ARIN6 NUMBER, 14X,17HP1V0T ANGLE RATI0//14X , 5HRI /RO , 22X , 1 4HBETA ! DE 
2GREES),19X,7HLAMBDA , 20X , IIHTHE TAP/BE TA//) 
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930 FORf'ATC 15X,G1D.S»2DX,G10.5,?0X»G10.5,213X,G1Q.5///) 

9*»C FORHATCSOX, 17HNUHERICAL RE S UL TS // 5X » 4 H ( R 1 ) , 8 X , 4H ( R 2 ) * 8X , 4H ( R 3 ) » 8 X , 
14H( R4) ,8X ,4H(R5 ) ,8X ,4H«R6 ) ,8X ,4H(R7 ) , 8X t 4H( R8) ,8X,4H (R9> ,8X ,5H< Rin 
2»,7X»5H(R11)//) 

950 FORKAK 5X,5HH1/H2,7X,4H1/H2»5X, 10HEPSIL0N/H?t3X,9HL0AD,WBAR,3X,9HU 
INIT-L0AD,2X, 10HFRICTI0N,F,4X,6HF/WBAR,5X»8HH2/H(CP),5X, 5HRtCP),5Xt 
210HTHETA(CP)/,4X,5HX(CP)/117X,4HBETA//} 

960 FORMAT! 11 (2X,G10.5) /> 

C 

STOP 

END 



o o o 
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TPF5(0) .XCART 

SUBROUTINE EUCL ID (NR »NA , R1 , BE TA , TKR AT ) 

..SUBPROGRAM PRODUCES COMMON BLOCK CONTAINING GEOMETRIC PARAMETERS.. 

COMMON/GEOM/ANGB,tHPR ,DR fDA 
C 

FACT=6. 2831853/360. S 2 PI RADIANS/360 DEGREES... 

ANGB=BETA*FACT BBETA IN RADIANS.. 

THPR:THRAT*AN6B aTHETAPiPiVOT ANGLE) IN RADIANS.. 

DR:( l.-RI)/FLOAT(NR) aRAOIAL INCREMENT.. 

DA=ANGB/FLOAT(NA) BANGLE INCREMENT IN RADIANS.. 

C 

RETURN 
END 



on no oo nooooo 
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TPFt(O) .XBEGN2 

SUBROUTINE BEGIN2 (E ,HALT,HMIN ,HRA TIO ,ERATIO |RI ,$ I 
C 

common/geom/btr»thpr,dr,da 

THETAPrTHPR 

BETA=BTR 

RAD90=1 .5707963 S90 DEGREES IN RADIANS.. 


..CALCULATE MINIMUM FILM RATIO AND BEARING NUMBER.. 


..REGION NUMBER 1... .PIVOT ANGLE = THETAP .LT. 0 ... 
10 IF( THETAP. GE. 0.0) GO TO 2D 

EPS = ERATIO/( l.-E RATI 0»SIN( THE TAP-BET A) » 
HMIN=EPS/ERATIO 

HRA7I0=( 1 .♦EPSOR I ♦SIN* THETAP > )/HMIN 
IF(HMIN.LT.HALT) RETURN 7 
GO TO IDO 


..REGION NUMBER 2... .THETAP ,GE. ZERO AND .LE. BETA ... 

20 IF(THETAP.GT,R£TA) GO TO 30 

EPS = ERATIO/ ( 1 ,-ERATIO«SlN(THETAp-BETA ) ) 

HMIN=EPS/ERA TIO 

HRATIO = ( 1 .♦EPS^SIN* THETAP) )/HMIN 
IFCHMIN.lt. HALT) RETURN 7 
GO TO 100 

..REGION NUMBER 3. ...THETAP .GT. BETA AND .LT. 90 DEGREES ... 
30 IFC THETAP, GT.RAD9OG0 TO MD 

EPS=ERATIO/{ 1 ,-ERATlO*RI«SIN CTHETAP-BETA ) ) 
HMIN=EPS/ERATIO 

HRATIOrC 1 .•»EPS«S1N(THETAP) )/HMlN 
IFCHMIN.LT.HALT) RETURN 7 
GO TO 100 


RO 


REGION NUMBER < 4 .... THETAP .GT. 90 DEGREES ... 

EPS = ERATIO/C 1 .-£RATIO^RI*SIN*THETAP-BETA ) ) 
HMIN=EPS/ERATIO 


HRATIO-C 1.+EPS)/HMIN 


IFCHMIN.lt. HALT)RETURN 7 


C 

100 CONTINUE 
E=EPS ■ 
RETURN 


END 


C 


on on ooooo 
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♦ TPFf < D» .COEFF 

SUBROUTINE ARR A Y S ( N R , NA , EPS , K K , RI , HL MBD A ) 

..SUBPROGRAM CALCULATES VALUES OF THE NODAL COEFFICIENTS AND THE INITIAL 
..VALUES OF THE DEPENDENT VARIABLE FOR THE FIRST ITERATION 

LOGICAL DEBUG, OLDO 

COMMON/GEOM/ANGB,THPR,DR,DA 

COMM ON /BLKA/A( 15,21 )/PLKB/B < 1 5 , 2 1 ) VB LKC /C f 1 5 , 21 ) /BLKD/D ( 15,21 ) 
C0MM0N/fiLKE/E<15,21 )/BLKG/G (15,21 >/PLKH/H ( 15 , 21 ) /BLKC/C (15,21) 
COMMON/BLKR/R(15)/ELKTH/TH(21 ) 
COMMON/BLOGIC/DEBUG,TABOUT,OLDQ,VARGlf;D 


LASTR=NR*1 
LASTA = NA-»1 
CR1= .5/OR 
CR2=1./DR**2 
CA1= .5/DA 
CA2r i./DA**2 
GO TO 5 


ENTRY RARRAY(EPS,KK ,HLMBDA) 

5 THETA = 0.D 

DO 2 JA=1,LASTA 
RADrRI 

TH(JA)=THETA 

ANG: (THPR-THETA) BANGLE FROM PIVOT LINE IN RADIANS 

STRIG=EPS*SIN( ANG) 

CTRIG=EPS*COS( ANG) 

C 

DO 1 JR=1,LASTR 
R(JR)=RAD 
HRA=1 .♦RAD*STRIG 
C 

IF(KK.EQ.O) GO TO 6 
IF( .NOT. OLDO) GO TO 6 
IF( JA.EC.l.OR.JA.EQ.NAPl) GO TO 6 
IF( JR.EO. I.OR.JR.EQ.NRPI) GO TO 6 
C 

HRATI0=HRA/H( JR,JA) 

Q(,jR, JA )=0( JR, JA)*HRATI0**2 
GO TO 7 
• C 

6 0(JR, JA)rHRA*>C2 

7 H(JR, JA )::HRA 

F.RST:1 ./RAD-STRIG/HRA 
SCND=CTRI6/(RAD*HRA ) 

CA2RAD=CA2/RADX'*2 

C 

A ( JR, JA )=CR2+CR1*FRST 
B ( JR, JA)-=CR2-CR1*FRST • 

C (JR, JA)=CA2RAD*CAltSCND 
D(JR, JA)=CA2RA0-CA1*SCN0 
E(JR, JA) = 2.*(CR2’^CA2RAD) 
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G( JR, JA)=HLMBDA*CA1/HRA 
C 

IF( . NOT. DEBUG )G0 TO 4 

WRITE (6,3) JR, JA, THETA, RAD flBEGIN DEBUG SEARCH 

WRITE (6,3) JR, JA,STRIG,CTRIG 

WRITE (6,3)HRA,Q( JR, JA ), A( JR,JA ) , B ( JR , JA ) , C ( JR , JA ) ,0 ( JR , JA ) , 
1 E (JR, JA) ,G( JR,JA) 

C 

4 RAD=RAD*DR 

1 CONTINUE 
THETA=THETA+OA 

2 CONTINUE 

3 FORMAK ) aFORMAT-FREE WRITE STATEMENT FOR DEBUG 

RETURN 

END 



n o r» o o 
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♦TPFt(O) .SEIDEL 

SUBROUTINE 


RELAX ( NR ,N A (OMEGA , ITERMX , RESIDE) 




..GAUSS-SEIDEL ITERATION WITH CONVERGENCE WHEN LARGEST DIFFERENCE BETWEEN 
..SUCCESSIVE ITERATIONS FOR ANY MATRIX ELEMENT IS LESS THAN THE INITIALIZED 
..PARAMETER "RESIDL" 

C0HM0N/BLKA/A(ir(21 ) /BLKB /B Cl 5 , ?1 I / BL KC / C ( 1 5 , 2 1 ) /BLK D/D { 1 5 , 2 1 ) 
C0HMCN/[»LKE/E(15,21 ) /BLKG/G ( 1 5 , 21 I /PLKH /H ( 1 5 , 21 ) /3LK Q/C> C 1 5 , 2 1 ) • 

COMMON/ELKP/PNOPM(15 ,21 )/BLKR/R ( 1 5 ) /BLK T H/ T H < 2 1 ) 

COMMON/BLCGIC/DFBUG , TABOUT,OLDQ,VARGRD 


C 


LOGICAL DEBUG, TABOUT 
REAL OMEGA 


30 

35 

36 


<«0 


FORMATC ) SFORMAT FREE WRITE FOR DEBUG 

FORMATC 1H0,3RX,A3HARRAY OF NORMALIZED PRESSURE PNORM ( ROW , COL ) // 
140X ,23HRADII FROM INNER R AD lUS/ /22X , 1 IF 1 0 . 3 / / ) 

FORMATC IHO, 1X.6HTHE TAr,F10.4, 3X , 2HP= , 11 F 1C) . // ) 

60 TO ‘♦C 

ENTRY RRELAX 
DO 3DD KK = 1 (ITERMX 
STRERRrO.Q 


DO 100 11=2, NR 

STORQrOC II , JJ) B Q-VALUE AT LAST ITERATION. 

GSOT=G (II, JJ ) /SORT CS TORO ) 

CCOF=CCII, JJ)-GSQT 
DCOF=OCII, JJ)+GSQT 

OCII , JJI = C l.-0ME6A)*ST0RQ + 0MEGA«M A C II , J J ) *0 ( 1 1 + 1 , J J ) + 
B(II,JJ)*QCII“1,JJ)*CC0F*0CII,JJ+1 )♦ 
DCOFOOCII ,JJ-in/ECII ,JJ) 

ERROR=OCII ,JJ)-STORQ 


SO 


IF ( .NOT.DEBUGIGO TO 50 

WRITEC6,3C)II,JJ,ST0R0,ERR0R,STRERR 

TEST=ABSCERROR) 


DEBUG 

DEBUG 


SWITCH 

TRACE 


C 

c 


100 

200 

300 

<t00 


4<i0 


IFCTEST.LE .STRERRIGO TO 100 
STRERR=TEST 
ISTORE=II 
JSTORE=JJ 
OSTORE=OCII, JJI 
CONTINUE 
CONTINUE 

IFCSTRERR.lt. RESIODGO TO 400 
CONTINUE 
CONTINUE 

IFC .NOT .DEBUGIGO TO 440 

WRITEC6 ,3C)KK ,STRERR,OSTORE,iSTORE,JSTORE,RESIOL 

LASTR=NR+I 

LASTArNA+1 


aUPDATE LARGEST RESIDUAL 
SROW OF LAST UPDATE 
aCOLUMN OF LAST UPDATE 
aO-VALUE AT LARGEST RESIDUAL 
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IF(DEBUG)WRITEI6,30)KK,STRERR,QST0RC,IST0REf JSTORE SDATA AT END OF 

aiTERATIONS. . . 

C 

DO <160 JA = 1,LASTA 
DO 450 JR=1,LASTR 

RNORM(JR,JA>=S0RT{0(JR,JA)/H(JR,JAI**2)-l. 3N0RMALIZED PRESSURE 

8RELATIVE TO ATMOSPHERE 

450 CONTINUE 
460 CONTINUE 
C 

IFI .NOT .TAB0UT)G0 TO 500 aSKIP NORMALIZED 

a PRESSURE PRINTOUT 

WRITE{6,35)(R<K l,K = l,LASTR) 

C 

DO 470 L=1,LASTA 

URITE(6,36MTH(L (PNORH(M,L }fH=l«LASTR) ) 

470 CONTINUE 
500 CONTINUE 
C 

RETURN 

END 
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rPFS(O) .PADCTR 

SUBROUTINE TABULT ( NR ,NA ,R I ,EX IT t L AMBD A, HMIN ) 

C 

C.. EXECUTIVE FOR INTEGRAL CALCULATIONS TRYING THE SIMPSON 1/3 RULF 
C..FOR ANGULAR COORDINATES... 

C 


C 


c 


1 


c 


5 


c 


c 


10 




REAL LAMBDA 

COMMON/CEOM/ANGP, THPRtDRtDA 

C0MM0N/BLKXXX/XF(21,3) fiARRAY OF INTEGRATION RESULTS FROM SUBR. RSIMP 
DIMENSION ODDiai.EVENCBl rENDI*!) ,EXIT(5) 

NA2=NA/2-l , 

LASTArNA+1 
GO TO 1 


ENTRY RTAB(EXIT, LAMBDA, HMIN) 

CALL RSIMP(NR,NA,LAMBDA,HMIN) SRADIAL INTEGRATION SUBROUTINE 


DO 5 1=1 

END( I) =0.0 
E VEN (I )=n.o 


ODDI I) =0.0 
CONiTINUE 


a SET STORAGE VECTORS TO ZERO.. 


ENOm=END{ 1)4XF(1,1)*XFILASTA, 1) 
END(2) = END( 2 ) ♦XF ( 1 , 2 ) +XF C L A ST A , 2 ) 
END( 3) = END( 3 ) ♦ XF < L A ST A , 2 ) *S IN ( A NGB I 
END(‘*)=END(4 )*XF(1, 3)*XF(LASTA,3) 


DO 10 JA=1,NA2 aiNTEGRATION ON INTERIOR OROS. 

EVENd ) = EVEN(n^XF(2#JA,l) 

F VEN (2)=EVEN (2)+XF (2*JA,2) 

EVEN (3)=EVEN(3) + XF «2*JA ,2)*S1;NI0A*FL0AT<2*JA-1 ) ) 

EVEN (4 )=EVEN (4 >+XF »2*JA,3) 

ODDI 1)=0DD(1)+XF(2.>«‘JA>1 ,1) 

0DD(2»=0PD{2)*XF(2*JA+1,2) 

ODDI 3)=0DDI3»*XF|2*JA*1 ,2)*SINI0A*FL0AT(2«JA) ) 

ODDI 4) =0DD(4 ) + XF 12+JA+l, 3) 

CONTINUE 

EVEN(1)=EVENI1)>XF(NA,1) slast even ordinate 

EVEN I2»=EVENI2)+XFINA,2) 

EVEN (3) =EVENI 3)+XFI NA ,2)*SIN(0A*FL0ATINA-1) )* 

EVENI4)=EVENI4 )+XFINA,3) 


WW=DA*IENDI n + 4.*EVENI 1 )+2.*0DD( 1 ) )/3. 
RR=DA*IEND( 2)*4.*EVENI2)+2,*0DDI2) )/3. 
AA=DA«IENDI 3)+4.YEVENI3)*2.*0DDI3) )/3. 
FF=DA*IENDI4 )*4.«EVENI4)*2.*0DDI4) )/6. 


aLOAD INTEGRAL. 
aRAOIAL MOMENT. 
aANGLE moment. 
aFRICTION INTEGRAL. 


AA=( ASINI AA/RR) i/ANGB 
RR=RR/WW 

WUNIT=2.*WW/IANGB*I 1.-RI+42)) 


aANGULAR C.P. COORDINATE. 
8RADIAL C.P. COORDINATE. 
aUNIT LOAD. 


EXITIl )=RR 
EXIT(2)=AA 
EXITI3)=WW 


31 


EXIT(4 )=FF 
EXIT(5I=UUNIT 

RETURN 

END 



no nnonoooo noon 
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•.♦TPFSCO.SIMPSN 

SUBROUTINE RSIHP ( NR , N A ,L AMBDA ,HMIN ) 
..SIMPSON INTEGRATION METHOD ALONG PAD RADII. 


REAL LAMBDA 

C0HM0N/GE0M/AN6P ,THPR ,DR,DA 

COMMON/BLKP/P (15,21 ) /BLKR/R ( 1 5 ) /8LKH/H( 1 5 ,2 1 ) 
COMMON/BLKXXX/XF (21 ,3 » 

DIMENSION END(3) ,EVEN(3),0DD( 3) ,CC(5,3l 

DATA CC/-1. i*!. ,-3. , 3*0.0, 1. ,0.D,-1.,3*D.0,3.,-‘».,1./ 


..USAGE OF FUNCTION DEFINITIONS .. 

F .INTEGRAND FOR LOADS AND MOMENTS, 

H? RATIO OF FILM THICKNESS TO MINIMUM FILM THICKNESS, 
X INTEGRAND FOR FRICTION MOMENT CALCULATIONS .. 


DEFINE H2(JP,JA»=H(JR,JA»/HMIN 
DEFINE F(K , JR, JA )=P (JR, JA »*R( JR )**K 

DEFINE X(M, JR,JA )-(CC (1 ,M)*P( JR , JA-»2) +CC (2,M)*P( JR ,JA>1 >♦ 

1 CC( 3 ,M )*P( JR, JA) ♦CC(4,M)*P( JR, JA-1 )♦ 

2 CC(5,M)*P(JR,JA-2) ) *R ( JR ) *H2 ( JR , JA ) / ( 2 . *DA ) ♦ 

3 LAMBOA*R ( JR)**3/ (3.*H2( JR, JA n 



NENDrNR/2-1 
LASTR=NR*1 
-LASTA=NA+1 
M = 1 

DO 50 JA=1,LASTA 

IF(JA.EQ.LASTA) M=3 

c ■ 

DO 10 1=1,3 9 SET STORAGE VECTORS TO ZERO.. 

END(I)=0.0 
EVEN(I)=0.0 
ODOd ) = 0.0 
IQ CONTINUE 

C 

END( 1)=ENDCU+F( 1, 1,JA)+F( 1,LASTR, JA) 

END(2)=END(2 »+F(2,l,jA)*F(2,LASTP, JA) 

END( 3)=END(3)*X(M, 1, JA)+X(M,LASTR, JA ) 

C • • 

DO 40 JJ=1 ,NEND s' 

EVEN( 1)=EVEN( 1)4F( 1,2*JJ, JA) 

EVEN( 2)=EVEN( 2)+F(2,2*JJ, JA ) | ^ 

EVEN( 3)=EVEN( 3)+X(M,2*JJ, JA ) 

ODD( 1 )=ODD( 1 ) +F ( 1 ,2*JJ*l , JA ) 
0DD(2)=0DD(2)^F(2,2*JJ*l, JA) 

ODD (3 )=00D( 3) ♦X(M,2*JJ*1, JA ) 
i»0 CONTINUE 

EVENd )=EVEN(1)^F(1,NR,JA) aLAST EVEN ORDINATE 

EVEN(2)=EVEN(2)^F(2,NR,JA) 

EVEN(3)=EVEN(3)*X(M,NR,JA) 
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00 3C 11=1,3 

XF(JA,II )=OR*<ENOf II)«4.»EVEN(II)42.*000(1I) ) /3. 
30 CONTINUE 
C 

M=2 

5CJ CONTINUE 
C 

RETURN 

END 
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TPFSIO.ELT 

SAMPLE DATA INPUT AS READ BY SUBROUTINE ‘MAIN*.. 

♦CARD 1-FORMAT 3I 6 ♦ F8 . 2 , 2E 8 . 1 NR' -NUMBER RADIAL MESH INCREMENTS, 

NA -NUMBER ANGULAR MESH INCREMENTS, 
ITERMX-MAXIMUM ITERATIONS IN GAUSS- 

SEIDEL ROUTINE IF CONVERGENCE 
FAILS, 

OMEGA -RELAXATION FACTOR IN GAUSS- 
SEIDEL METHOD, 

HALT -MIN LUBRICANT FILM THICKNESS.. 
RESIDL-CONVERGENCE C R I TERION= TE S T ON 
MAXIMUM CHANGE IN Q-VARIABLE 
BETWEEN ITERATIONS.. 


NR 

8 


NA ITERMX OMEGA HALT RESIDE 

12 ROD 1.25 .lE-3 .lE-4 


♦CARD 2-FORMAT 4L6 DEBUG (.TRUE.= PRINTOUT OF COMPLETE ARRAY OF 

DEPENDENT Q-VARIABLE FROM 
6AUSS-SEIDEL ROUTINE AT EACH 
ITERATION) , 

TABOUT (.TRUE.: PRINTOUT OF PRESSURE ARRAY 
OVER PAD AREA AT 
CONVERGENCE), 

. OLDO (.TRUE. : USE CONVERGED VALUES OF Q- 

VARIABLE AS STARTING ESTIMATES 
FOR NEW CALCULATIONS WITH 
DIFFERENT INITIAL CONDITIONS). 

VARGRO (i,TRUE.: CREATES VARIABLE MESH GRID 
OVER SECTOR PAD FOR LIQUID 
FILM CALCULATIONS) . 


’BUG TAB* OLD* VAR'D 
FALSE FALSE TRUE TRUE 


♦CARD 3-FORMAT 5110 ..INDICES FOR PARAMETER ARRAYS IN NAMELIST 

INPUT ''VAR BLE" i,, 

NTHETA NRATIO NUHRI NUHBET NUMLKB 

1 20 1 1 1 


♦CARDS 4,5,,, .;DATA INPUT FOR NAMELIST "VARBLE" .,. 

...DESCRIPTION OF INPUT PARAMETER ARRAYS ... 

NUMBER OF PARAMETER PARAMETER DESCRIPTION 

DATA ARRAY 
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******** ******** *^**************** 

ntheta tratio pad pivot angle, thetap/beta 

NRATIO eratio pad slope factor .epsilon/hmin 

NUMRI VRl PAD INNER RADII, RI 

NUMBET vbeta pad angle dimension, beta 

NUMLMB VLMBDA BEARING NUMBER, LAMBDA 


...INPUT FOLLOWS ... 

SVARBLE TRATI0=1.0,VPI:.5,VBETA=45.,VLMBDA=50., 

ERATIO=.5,l.,l.S,2.,2.5,3.,3.5,4.,4.5,S.,5.S,6.,6.5,7..7.5f8., 
8.5,9. ,9.5,10., ' 


SEND 
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FSTBP2 

Flat Sector Pad Thrust Bearing E^rogram Number 2 — liquid lubricant. 


*TPFS( G » .M AIN2 
C 


C**HAIN EXECUTIVE FOR GAUSS-SEIDEL ITERATION, . LIQUID LUBRICANT., 
C 

LOGICAL DEBUG, TABOUT,OLDQ,VARGRD 

c 

COMMON/BLOGIC/DFBUG ,T ABOUT, OLOQ fVARGRD 

COMMON/GEOM/ANGB ,THPR ,0R,DA 

DIHENSI ON EXIT ( 9) ,TRATI0( 15 ) ,ERATIO (26) 

DIMENSION VRK 1C) ,VBETA( 10) ,VLMBOA( 10) 

C 

NAMELIST/VARBLE/TRATIO,ERATIO,VRI,VBETA,VLMBDA 

C 


1 

2 

3 

4 


5 

C 

c 


c 


c 


c 


c 

c 


c 

c 


8 

ID 


C 


FORMAK 3I6,F8.2,2E8.1 ) 

FORMAK 4L6) 

F0RMAT(5110) 

FORM AT( 2X , 12HN0 . OF RCWS=,I3,16H NO. OF C0LS = ,I3,26H MAX NO. 
1 OF ITERATIONS: ,14, 25H RELAXATION PARAMETER: , F4 .?, 19H RESIDU 
2AL ERROR:, Fia.7/2X,38HSMALLEST ALLOWED FILM THICKNESS HALT r, 
3610,4//) 

F0RMAT(2X,39HTHF MINIMUM FILM THICKNESS IS LESS Th AN ,G10 . 4 / / ) 


READ <5, 1) NR, N A, ITER MX, omega ,H alt, RESIDE 
READ (5, 2)DEBUG,TABOUT,OLOO, VARGRD 
READ <5, 3)NTHETA,NRAT10,NUMRI,NUMBET 
READ (5, VARBLE) 


WRITE! 6,4 )NP,NA ,ITERMX,OMEGA,RE 

WRITE(6,900) 

IPRINT:0 

00 500 NRI:l,NUMRI 
RIrVRKNRI) 

DO 400 NBETA:1 ,NUMBET 
BETA:VBETA(NBETA ) 


IDL,HALT 8NUMERICAL SOLUTION 
apRINTOUT... 
a HEADING PRINTOUT., 
a INITIALISE PRINT SWITCH... 

a INCREMENT INNER RADIUS VALUES.. 


a INCREMENT BETA VALUES , .DEGREES. 


DO 200 NANGLE:1,NTHETA a INCREMENT THETAP/BETA RATIO 

XXTHTP:TRATI0(NAN6LE) 

CALL EUCLIDtNR,NA,RI,BETA,XXTHTP,VABGRD) aCALCULATES GEOMETRIC 

PARAMETERS. . . 

KOUNT:0 

IF(IPRINT,GT,0)G0 TO 8 
WR1TE(6,910) 

GOi TO ID 

WRITE(6,915) 8PRINT START TOP OF PAGE.. 

«RITE(6,920) 

WRITE(6,930)PI ,BETA,XXTHTP 
HRITE(6,940) 

WRITE(6,950) 

DO IOC NFILMri ,NRAT10 8 INCREMENT EPS/H2 RATIO 

ERAT:ERATI0(NFILM) 


o o n o noo o n non o n nn 
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CALL BEGIN^ICPSfHALTtHHINfHRAT t£RAT,RItS12) 
GO TO 15 


3FILM THICKNESS 
9 CALCULATIONS.. 


12 WRITE(6,5>HMIN 9 ERROR EXIT.. 

60 TO 100 

15 IFIKOUNT.GT.OIGO TO 21 

..THREE SUBROUTINES TO SET UP AND SOLVE THE FINITE DIFFERENCE 
..REYNOLDS EQUATION, AND RETURN INTEGRALS OF PRESSURE OVER PAD AREA.. 
CALL ARRAYS! NR,NA,EPS,KOUNT,RI ,HNINI 
CALL RELAX(NR,NA ,ONEGA ,ITERMX,RESIDL, VAKGRDI 
CALL TABULT(NR,NA,HMIN,RI, EXIT, VARGRD.EPS) 

60 TO 22 


21 


CALL RARRAY<EPS,KOUNT,HMIN ) 9ENTRY TO SUBROUTINE ARRAYS. 

CALL RRELAX 90ITT0 FOR GAUSS-SEIDEL ROUTINE.. 

CALL RTAB(HMIN,EXIT) 9ENTRY TO INTEGRAL CALCS.. 


22 


Rl = 

ERAT 

R2 = 

HRAT 

R3 = 

EXIT!7) 

R4 = 

EXIT!3»/EXIT 

R5 = 

EXIT!4) 

R6 = 

EXIT !55 


9RATI0 EPS/HMIN «. 
9RATI0 H1/H2 .. 
aUNIT LOAD .. 
9FRICTI0N/L0AD .. 
9R0 SIDE leakage .. 
9RI SIDE LEAKAGE .. 


♦♦R7 IS TRAILING EDGE LEAKAGE,,, 
**R8 IS FLOW INTO LEADING EDGE ... 
R7|: EXIT(9) 

R8= EXIT(8) 


R9=IEXIT(2)-PI I/(1.-R1) 

R10= EXIT(6) 9THETA C.P. .. 

Rll= EXIT!2M'SIN( ANGB*R10-THPR) 


WRITE! 6,960 )R1,R2»R3,R4 tR5,R6,R7,R8,R9, RIQ,R 11 

c ■ ■ ■ • » . 

K0UNT=K0UNT+1 

100 CONTINUE • • 

, ■ . IPRINT = IPR1NT^1 

200 CONTINUE 
C 

400 CONTINUE 
500 CONTINUE 
C 

900 FORMAT! 1H1//40X,43HNUMERICAL RESULTS - INCOMPRESSIBLE ANALYSIS/ 
140X , 28HLEWIS RESEARCH CEN TER ! N AS A ) / / / / > 

910 FORMAT! BOX, 17HSYSTEM PARAMETERS // 1 5X , 4H ! P 1 » ,25X,4H!P2) ,28X ,4H!P3 ) , 
127X,4H!P4)//> 

915 FORMAT! IHl / /50X , 17HSYSTEM PAR AMETERS/ /1 5X, 4H <P1 ) , 25X , 4H !P2 » , 28X , 
14H(P3),27X*4HCP4)//) 

920 FORMAT! lOX, 17HINNER/0UTER RADII, 14X,17HBEARING PAD ANGLE ♦ 14X, 17HPI 
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JVOT angle RATI0//14X,5HRI/R0, 22X, 14HBET.A ( DEGREE S ), 1 9X , 1 IHTHETAP/B 
2ETA//) 

930 F0RNAT( 15X,Gin.5,20X,G10.5,20X,61D.?/// ) 

940 FORMAK BOX, 17HNUMER1CAL RESULTS//BX,4H(R1), 8X,4H ( R2 ) , 8X , 4H ( R3 ) , 8 X, 
14H( P4) , 8X ,4H tRB ) ,8X ,4H(R6 » ,8X ,4H<R7) ,8X ,4H»R8) ,8X,4H IR9 ) ,8X,5H(P 10 
2) ,7X ,5H (Rll )//) 

950 FORNAK 1 X , 1 THEP S I LON /H2 , 5X , 5HH 1 /H2 , 5X , 9HUNI T-L 0 AD , 4 X , 9HFRI C T I ON/ , 
14X,7HR0 SIDF. ,4X,7HRI SIDE , 6 X , 4HE0GE , 8 X , 6HI N S IDE , 3X , 1 IH ( R ( CP ) -R1 ) / , 
22X, ICHTHETA (CP )/ ,4X,5HX(CP) /44X .4HL0AD, 4X, 7HLEAKAGE, 4X, 7HLEAKAGE , 
~J5X, 7HLEAKAGE.7X ,4HFL0W,7X , 6H ( 1 -R I ) , 1 OX , 4HBE T A // ) 

960 FOR^'AK 11 C2X,Gin.5) /) 

970 FORMATt 3(G15.5M 
C 

STOP 

END 


U O O U U (J 
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TPFS(0).XCART 

SUBROUTINE EUCLID (NR ,NA ,R I , BE TA , THR A T , V ARGRD ) 

..SUBPROGRAM PRODUCES COMMON BLOCK CONTAINING GEOMETRIC PARAMETERS.. 
LOGICAL VARGRD 

COMMON/GEOM/ANGP ,THPR ,DR i D A /BLK R/R ( 2 3 ) , DELR ( 23 > 

COMMON/BLKTH/TH C25 ) , DELTH ( 2 5 » /BLK S I N / XS I N ( 2 5 ) , TS IN I 2 5 ) » TCOS ( 25 ) 
COMMON/INDEX/LN(23) , L S ( 2 3 ) , LE ( 2 5 ) , LW ( 25 ) 

INTEGER HALERfHALFA 

FACT=6. 2831853/360. 

ANGB-BETA*FACT 
THPR=THRAT«ANGB 
DR= ( l.-RI )/FLOAT(NR) 

OA=ANGB/FLOAT(NA) 


INTE6R = 1 
IFCVARGRD ) INTEGRA? 

L ASTR:NR-*-INTE6R 
LASTAtna+INTEGR 
HALFR=NR/2+1 
HALFA=NA/2+1 

IF(VARGRD)GO TO 1 
DELTH(l) = DA 

DELTH(LASTA> = DA 
OELR(l) .= DR 
DELRCLASTR) = DP 
60 TO 2 


1 

DR*» 

r 

.25*DR 


DA4 

— 

.25«DA 


DR2 

- 

2.«0R4 


DA2 

z, 

2.*DA4 


DELP( 1 ) 

z 

DR4 


DELTH( 1) 

z 

DA4 


DELR(LASTR) 

z 

DP4 


DELTH(LASTA) 

'Z. 

DA4 


C 

2 R(1)=RI 


R(LASTR)=1 

TH(1)=0.0 

TH(L ASTA )=ANGB 

XSIN ( 1 ) =0.0 

XSIN (LASTA ) rSIN( ANGB) 

TSINd ) r SIN(THPR) 

TSIN (LASTA )=SIN( THPR-ANGB) 

TCOS ( 1 ) = COS( THPR J 

TCOS (LASTA ) =COS( THPR-ANGB) 

C.. INDICES FOR VARIABLE MESH DIFFERENCE EQUATIONS. .. 
KPLUSn 
MINUSrLASTR 
C 

LE(l):i 

LW(LASTA )=LASTA 


a 2 PI RADIANS/360 DEGREES... 
dBETA IN RADIANS.. 
aTHETAP(PIVOT ANGLE) IN RADIANS.. 
3RADIAL INCREMENT.. 

BANGLE INCREMENT IN RADIANS.. 


40 


C 


C 

10 


c 

20 


V, C 

30 


50 

C 

100 

c 


c 



c 

110 


c 

120 


C 

130 


LS(1> = 1 

LN(LASTR )=LASTR 

DO 100 jritHALFR 
K=1 

DDR=DR 

IF(.N0T.VARGRD)60 TO 30 
IFCJ.GT.2) GO TO 30 
60 TO { 10,20 » , J 

K=3 

DDRrDR** 

NCWri 
60 TO 30 

Kr2 
N0U = 2 
DDR = DR2 

DO 50 L=1,K 

IFIJ.FQ. 3 )N0W=2 
KPLUS=KPLUS+1 
MINUSrLASTR+l-KPLUS 
LSIKPLUS ) =KPLUS-NOW 
LNCKPLUS ) =KPLUS*1 
LSCMINUS ) =MINUS-1 
LNCMINUS ) =MINUS+NOW 
DELR(KPLUS) : DOR 
DELR(MIRUS» = DDR 

R(KPLOS» = R(KPLUS-l>+DELR»KPLUS-n 
RCMINUSI = R(MINUS+1>-DELR(MINUS+1) 
N0W = 1 
CONTINUE 

CONTINUE 


KPLUS=1 

MINUS=LASTA 

DO. 200 J=1,HALFA 
K=1 

D0A=0A 

IFr.N0T..VAR6RD)G0 TO 130 
IF«J,GT.2) GO TO 1.30 
60 TO (110, 120), J 

K=3 

0DA=0A4 
Nown 
GO TO 130 

K=2 

00A=DA2 

N0Wr2 

DO 150 L=1,K 
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IF( J.E:Q.3)N0W:2 
KPLUS=KPLUS+1 ‘ 

MINUS=LASTA+1“KPLUS 
LW(KPLUS)rKPLUS-NOW 
LE(KPLUS)=KPLUS+1 
LE(HINUS)=MINUS-»NOW 
LW(MINUS)=MINUS-1 
DELTH(KPLUS) = DDA 
DELTH(PINUS) = DDA 

TH(KPLUS) = TH(KPLUS-1 ) ♦DELTH(KPLUS-1 » 
TH(MINUS) = TH(MINUS+1)-0ELTH IMINUS+1) 
TEMPI - TH(KPLUS) 

TEMP2\t: TH( MINUS) 

X SIN ( KPLUS ) rs IN (TEMPI ) 

XSIN(MINUS)= SIN(TEMP2) 

TSINiKPLUS )=S IN (THPR -TEMPI ) 

TSIN(MINUS)= SIN( THPR-TEMP2) 
TC0S(KPLUS)=C0S(THPR-TEMP1I 
TCOS( MINUS }=C OS CTHPR-TEMP2} 

NOW = l 

150 CONTINUE 

C 

200 CONTINUE 
C 

RETURN 

END 


r> o o o r> o 
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♦TPFf (0) .XBEGN2 

SUBROUTINE BEG I N2 (E » H AL T , HMIN ,HRA TIO , ER ATIO t RI 1 5 I 
C 

C0MMCN/GECM/AN6B,THPR.0RiDA 

THETAPrTHPR 

beta =ANGB 

RAD90=1 .5707963 390 DEGREES TN RADIANS.. 


..CALCULATE MINIMUM FILM RATIO AND BEARING NUMBER.. 


..REGION NUMBER 1... .PIVOT ANGLE = THETAP .LT. 0 ... 

10 IF( THETAP. GE. 0.0) GO TO 20 

EPS=ERATIO/ ( 1 =-ERATIO*SIN( THETAP-BETA ) I 
HHINrEPS/ERATIO 

HRATIO: < 1.4EPS^RI«SIN( THETAP ) ) /HMIN 
IF(HMIN.LT.HALT) RETURN 7 
GO TO 100 

..REGION NUMBER 2... .THETAP .GE. ZERO AND .LE . BETA ... 

20 IF( THETAP. GT. BETA ) GO TO 30 

EPS=ERATIO/( 1 .-ERATIO»SIN( THETAP-BETA ) ) 
HMIN=EPS/ERATIO 

HRATIO=( 1 .♦EPS*S IN (THETAP ) )/HMlN 
if(hmin.lt.halt) return 7 

GO TO 100 

..REGION NUMBER 3. ...THETAP .GT. BETA AND .LT,. 90 DEGREES 
3U IF( THETAP. GT.RADPC) GO TO 40 

EPS = ERATIO/ ( 1 .-ERATIO«RI»SIN (THETAP-BETA ) ) 
hmin=eps/eratio 

HRATI0=( l.+EPS^SIN(THETAP) )/HMIN 
IF(HMIN.LT.HALT) RETURN 7 
GO TO 100 

..REGION number 4.... THETAP .GT. 90 DEGREES ... 

40 EPS=ERATIO/(l.-ERATIO*Rl*SlN(THETAP-BETAn 

HMIN=EPS/ERATI0 
HRATI0 = ( l.+EPS)/HM),'N 
IF(HMIN.LT.HALT)RC:TURN 7 


100 


CONTINUE 

E-EPS 

RETURN 


END 


noon 
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TPFS(OJ .COEFF 

SUBROUTINE ARP A YS ( N R , NA , EPS , K K , R I , HHI N ) 

..SUBPROGRAM CALCULATES VALUES OF THE NODAL COEFFICIENTS AND THE INITIAL 
..VALUES OF THE DEPENDENT VARIABLE FOR THE FIRST ITERATION . L lOUID LURE 

LOGICAL OLDO,VARGRD 

COM MON/ ELK A /A ( 2 3 , 25 ) /BLKB /B ( 2 3 , 25 ) /PL KC /C ( 2 3 , 25 ) /BLK D/D (23,25) 
COMMON/BLKE/E (23 ,25 ) /BLKF /F ( 2 3 , 25 ) /« L KH /H < 2 3 , 25 ) /BLK 0/Q ( 2 3 , 25 ) 
COMMON/BLKR/R (23 ) , DELR ( 23 ) /BL OG IC /DE b , T A B , OLDO , V A RGRD 
COMMON/ELKTH/TH (25) ,DELTH ( 2 5 ) /BLK SIN/ XS IN ( 25 ) , TS IN ( 25 ) , TCOS ( 2 5 ) 

C 
C 

INTEGRri 

IF( VARGRD )INTF6P = 9 
LASTA=NA*INTEGR 
LASTPrNRflNTEGR 
GO TO 5 
C 
C 

ENTRY RARRAY(EPS,KK ,HMIN) 

5 DO 2 JA=1 ,LASTA 
DA=DELTH(JA) 

CA1=.5/DA 
CA2=1 ./DA«»2 
STRIGrEPS^TSIN ( JA ) 

CTRI6=EPS*TC0S ( JA > 

C 

DO 1 JR=1,LASTR 
RAD=R(JR) 

DR=DELR(JR) 

CR1=.5/DR 
CR2:i./DR*«2 
RAD2=RAD**2 
HRA=1 .♦RAD^STRIG 
C 

IF(KK.EO.O) GO TO 6 
IF( .NOT.OLDQ) GO TO 6 
IF( JA.EO.l.OR.JA.EO.LASTA) GO TO 6 
IF( JR.EO. l.OR.JR.EO.LASTR) GO TO 6 
C 

GO TO 7 

C .• ^ • 

•6 Q(JR,JA) = 0.0 - 

7 H( JR, JA)::HRA 

RH = RADOHRA 

FRST = CRl-M^. + HRA-S.I/RH 
SCND = 3.*CTRIG*CA1/RH 
CA2RAD z CA2/RAD2 

C' ■ ' ■ ^ 

ACJR, JA)zcR2+FRST 
B(JR, JA)ZCR2-FRST 
C (JR, JA) ZCA2RAD-SCN0 
D( JR, JA) ZCA2RA0 + SCN0 
E ( JR,JAIz2.<c( CR24CA2RA0) 

F ( JR, JA)zCTRIG*RAO*HMIN**2/HRA**3 



CONTINUE 

CONTINUE 

return 



on on o on on nnono 
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*TPF?< 0) .SEIDEL 

SUBROUTINE REL AX ( NR . N A , 0ME6 A , ITERMX , RES I DL , V AR6R0 ) 


. .GAUSS-SEIOEL ITERATION WITH CONVERGENCE WHEN LARGEST DIFFERENCE BETWEEN 
..SUCCESSIVE ITERATIONS FOR ANY MATRIX ELEMENT IS LESS THAN THE INITIALIZED 
..PARAMETER "RESIDL" 

C0MM0N/BLKA/A(2 3.2S)/BLKB/B (2 3,25»/BLKC/C(2 3,25 )/BLKD/D (23,25 » 
COHMON/BLKE/E(23,25 )/ PLKF/F ( 2 3 , 25 ) /PL KH /H ( 2 3 , 25 )/3LKQ/Q (23,25 ) . 
COMMON/IN DEX/LN( 23) ,LS( 23 » ,LE (25) ,LW( 25 » 


LOGICAL VARGRD 
REAL OMEGA 


INTEGRri 

IF( VARGRD)INTE6R=9 

LASTR=NR+INTEGR-1 

LASTA=NA*1NTEGR-1 

ORIGINAL PAGE IS 
OB POOR QUALITY 


60 TO 4D 


ENTRY RRELAX 
40 DO 3 00 KK=: , TTERMX 
STRERR=0.0 
DO 200 JJ=2,LASTA 

IE=LE(JJ) a ANGULAR INDICES... 

IW=LW(JJ) 


DO 100 II=2,LASTR 

IN=LN(II) a RADIAL INDICES... 

ISrLSdl) 

STOROzO (II, JJ) a 0-VALUE AT LAST ITERATION... 

FIR ST=( 1 .-OMEGA )*STORO 
RADIALrA (II,JJHfO(lN,JJ)-»B(II ,JJ)»«0( IS,JJ) 

ANGULRrC (II, JJ)«0( II ,IE)*D (II ,JJ )<'OtII ,IW) 

Q( II, JJ) =FIRST+OMEGA*(RADIAL+ANGULR+F(II , JJ) )/E(II , JJ> 
ERR0R=0(II , JJ)-STORO 

50 TEST:A8S(ERR0R) 


IF(TEST.LE.STRERR)G0 to 100 

STRERRrTEST 

ISTORE=II 

JSTORE=JJ 

OSTORE-Q(II,JJ). 

100 CONTINUE 

200 CONTINUE 

IF(STRERP.LT.RES1DL)G0 TO 400 
300 CONTINUE 
400 CONTINUE 


a'uPDATE LARGEST RESIDUAL 
dROW OF LAST UPDATE 
aCOLUMN OF LAST UPDATE 
aO-VALUE AT LARGEST RESIDUAL 


RETURN 

END 


ono o ooon 
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♦TPFSI 0) .PADCTR 

SUBROUTINE TABUL T ( N R , NA , HMI N , RI ♦ EXI T , V AR GRD , EPS J 

..EXECUTIVE FOR INTEGRAL CALCULATIONS TRYING THE SIMPSON 1/3 RULE 
..FOR angular COORDINATES... 

C0MM0N/GE0M/AN6B , THPR,DR,DA 

COMMON/BLKXXX/XF ( 25,6) SINTEGRATION RESULTS FROM SUBR. RSTMP 
COMMON/BLKTH/TH (25) ,DTH ( 26)/BLKSIN/XSIN (25 ) ,TSIN (25) ,TC0S(25 ) 
COMMON/INDEX/LN (23) , LS ( 2 3 ) , LE ( 2 5 ) , L W ( 25 ) 

DIMENSION EXIT(9) ,SUHMA(7) 

LOGICAL VARGRO 

..FUNCTION DEFINITION FOR SIMPSON INTEGRATION RULE.. 

' DEFINE AREA ( J , K ) =DT H ( J* 1 ) ♦ ( XF ( J , K ) + *• .*XF ( J+1 ,K ) +XF ( J + 2 , K ) ) 
DEFINE TRI6(J) = DTH ( J* 1 ) ♦ ( XF ( J , 2 ) OXSIN ( J ) + 

1 <4.*XF( J+1,2)«XSIN( J+1) *XF( J+2,2)*XSIN( J + 2) ) 

C 

INTE6R=1 

IF(VAR6RD)INTEGR=9 

LASTArNA+INTEGR STRAILIMG EDGE NODE .. 

LASTR^NR+INTEGR 80UTER RADIAL NODE .. 

IQUI T: (LASTA-1 )/2 3 ANGULAR MIDPOINT .. 

FAC1=1.-RI**2 
FAC2= ( 1 .-RI«=C3)>CEPS 
C 

GO TO 1 
C 

ENTRY RTA8(HMIN,EXIT) 

1 CALL RSIMP(LASTR,LASTA,HMIN) 3RADIAL INTEGRATION SUBROUTINE 

FAC3= .5«FAC 1/HMIN 
FAC'4 = FAC2/ (3.«HMIN) 

C 

DO 2 1-1,7 3 SET STORAGE VECTORS TO ZERO.. 

SUMMA(I)=O.D 

2 CONTINUE 
C 

Cs. INDEX LIST FOR STORAGE VECTORS AND FOR RETURN VECTOR “EXIT",,, 


C 

(1 ) 

CALCULATIONS 

FOR 

TOTAL LOAD.. 



C 

( 2 ) 

H 


RADIAL MOMENT AND *C.P. 

COORDINATE. . 

c 

(3 ) 

M 


PAD FRICTION* . 



c 

(*») 


•1 

LEAKAGE FROM OUTER PAD 

ARC 

( R = 1 ) . . 

c . 

(5) 


u 

" ' •• INNER •• 

II 

(R=RI) * . 

0 

(6) 

• •• 


ANGULAR MOMENT AND C.P 

. COORDINATE .. 

c 

(7) 

H 


UNIT LOAD «. 



c 

(8 ) 


II 

FLOW INTO LEADING EDGE 

• • 


c 

(9) 


M 

TRAILING EDGE LEAKAGE 

• • 


c 

DO ID ISUK=1,IQUIT 


3SUM OVER PAIRS OF INTERl?ALSV« 


JFRD = 

2*ISUM-l 


30DD NODES.. 




DO 9 K 

= 1,6 


3TEMP0RARY STORAGE . . 




SUMMA (K)=SUMMA(K)-»AREA( JFRD,K) 

9 CONTINUE 
SUMHA(7)=SUMMA(7)+TRIG( JFRD) 

10 CONTINUE 



DO 12 KKK=1,2 


3EXIT VECTORS 


LL=*»*(KKK-1 ) 

HM=LL*1 

NN:MM^2 

DO 11 ILK=MM,NN 

ILK0rILK-(KKK/2 ) 

EXIT(ILKO) - SUMMAC ILK) /3. 

CONTINUE 

CONTINUE 

EXI T <6) = ( ASIN(EXIT( 6) /EXIT( 2) ) ) / ANCIB 
EXIT (2) -EXIT{2)/EXIT(1) 

EXIT(3)=EXIT(3)/6. 

EXIT <7) =2.*EXIT ( 1 )/ < ANGB+FACl > 
EXIT(8)=FAC3 + FAC4*XSIN(l)-XF(li‘») 

EXIT (9) =FAC3+FAC4*XSIN(LASTA)-XF<LASTA,4 ) 

RETURN 

END 



♦TPFS(O) 



C 


C 

C 

10 

c 


i{ 

!l 


15 
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.SIMPSN 

SUBROUriNE RSIMP(LASTR,LASTA,HMIN) 

SIMPSON INTEGRATION METHOD ALONG PAD RADII, LIQUID CASE.. 


COHMON/BLKO/P(2?, 25 )/BLKH/H<2 3, 25)/BLKR/R(2 3 ),DR(23I 
COMMON/ELKXXX/XF ( 25 ,6)/BLKTH/TH (25) ,DTH (25) 

C0MM0N/INDEX/LN(23) ,LS(23),LE(25) ,LW(25) 

DIMENSION SUMMA(4) ,CC(5,3) 

DAT A CC/”1. ,4. ,“3. , 3^0.0, 1. ,0.C,**1. ,3’^Q.O, 3. ,1 ./ 

USAGE OF FUNCTION DEFINITIONS .. 

F INTEGRAND FOR LOAD AND MOMENTS, 

H2 RATIO OF FILM THICKNESS TO MINIMUM FILM THICKNESS, 

Y PARTIAL DERIVATIVE OF PRESSURE WITH RESPECT TO THETA, 

X INTEGRAND FOR FRICTION MOMENT CALCULATIONS, 

7 '• •' CALCULATION OF LEADING EDGE C TRAILING EDGE 

FLUID LEAKAGE .. 

DEFINE H2(JR,JA)=H( JR,dA)/HMlN 
DEFINE F(K, JR,JA )=P(JR,JA )*R(JR )«*K 

DEFINE Y( M, JR, JA , JE , JW):: (CC ( 1 ,M )«P( JR ,JE-»1 ) ♦CC<2 ,M )*P( JR, JE) + 

1 CC( 3 ,M ) *P( JR , JA ) +00(4 ,M ) ♦P ( JR , JW ) ♦ 

2 CC(5 ,M)*P( JR , JW-i ) )/ (2.*DTH( JA ) ) 

DEFINE X{ M, JR, JW, JE , JW) = 3.*Y(M, JR , JA , JE , JW )*R(JR )«H2 (JR , JA ) ♦ 

1 R(JR)**3/H?( JR,JA) 

DEF I NE Z ( M , JR , JA , JE , JW ) : Y ( M , JR , J A , JE , JW ) *H2 ( JR , J A ) 3/R ( JR ) 


DEFINE AREAl (K , JR, JA )z DR ( JR + 1 ) ♦ ( F (K , JR , JA > *4 . ♦F ( K , JR + 1 , JA ) + 
1 F (K , JR + 2, JA ) ) 

DEFINE AREA2(M, JR, JA, JE, JW)=DR( JR*1)*(X( M, JR,JA, JE, JW) ♦ 

1 .‘(.♦X{M,JR + 1,JA,JE, JW) + 

2 , X(M,JR+2,JA,JE, JW) ) 

DEFINE AREA3CM, JR, JA , JE, JW) =DR( JR+1 )*(2( M, JR ,JA , JE , JW)+ 

1 4.*ZCM , JR^l ,JA, JE, JW) ♦ 

2 Z(M, JR-»2,JA, JE, JW) ) 
IQUIT=(LASTR-1 )/2 

M = 1 

DO 50 JA=1 ,LASTA 

•IF(JA.EQ .LASTA ) M=3 

00 10 1=1,4 a SET STORAGE VECTORS TO ZERO.. 

SUMHA(I):0,0 


CONTINUE 

JE=LE ( JA ) 

JW=LW(JA) 

DO 20 XSUM=i,ICUIT aSUM OVER PAIRS OF INTERVALS*. 

JR=2*ISUM-1 aOOD nodes.. 

DO 15 K=l,2 aSTORE RESULTS IN SUMMA VECTORS.. 

SUMMA (K )=SUMMA(K)+AREA1CK, JR,JA) 

CONTINUE . 
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SUKMA( 3) rSUHMA ( 3 ) ♦AREA2( 3 , JR » J A , JE , JW » 

1F(M,EQ. 1 .OR.H.E0.3)SUMMA(4)=SUMWA ( 4 ) ♦ ARE A3 ( M f JR » JA i JE i JW ) 
20 CONTINUE 

C 

DO 40 11=1,4 

XF(JA,II)=SUMMA(II)/3, 

40 CONTINUE 
C 

1 

1 
C 

c 

M=2 

50 CONTINUE 
C 
C 

RETURN 
END 


XF(JA,5) = .5*(4.'«‘P<LASTR-1, JA ) - 

P(LASTR-2, JA ) )*H2 (LAS TR,JA)*» 3/DR ILASTR-n 
XF( JA,6I = .5«P( 1) ♦( -P( 3, JA) ♦ 

4.*P(2,JA) )»H2(lfJA)*<‘3/0R(l) 
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KTPF*(Q),ELT 

SAMPLE DATA INPUT AS READ BY SUBROUTINE ‘MAIN*.. 

♦ CARD 1-FORMAT 216 ,F6 .2 ,2E8 . 1 NR -NUMBER RADIAL MESH INCREMENTS , 

NA -NUMBER ANGULAR MESH INCREMENTS, 
ITERMX-MAXIMUM ITERATIONS IN GAUSS- 

SEIDEL ROUTINE IF CONVERGENCE 
FAILS, 

OMEGA -relaxation FACTOR In GAUSS- 
SEIDEL METHOD, 

HALT -MIN LUBRICANT FILM THICKNESS.. 
RESIDL-CONVERGENCE CRITERION=TEST ON 
MAXIMUM CHANGE IN Q-VARIABLE 
BETWEEN iterations.. 

NR NA ITERMX OMEGA HALT RESIDL 

IR 16 t!00 1.25 .lE-3 .IE-** 


♦CARD 2-FORMAT RLB DEBUG (.TRUE.= PRINTOUT OF COMPLETE ARRAY OF 

DEPENDENT Q-VARIABLF FROM 
GAUSS-SEIDEL ROUTINE AT EACH 
ITERATION), 

TABOUT (.TRUE.r PRINTOUT OF PRESSURE ARRAY 
OVER PAD AREA AT 
CONVERGENCE) , 

OLDO (.TRUE. : USE CONVERGED VALUES OF 0- 

VARIABLE AS STARTING ESTIMATES 
FOR NEW CALCULATIONS WITH 
different INITIAL CONDITIONS). 

VARGRO (.TRUE.r CREATES VARIABLE MESH GRID 
OVER SECTOR PAD FOR LIQUID 
FILM CALCULATIONS). 


•BUG TAB* OLD* VAR*0 
FALSE FALSE TRUE TRUE 


♦CARD 3-FORMAT 5110 ..INDICES F0R PARAMETER ARRAYS IN NAMELIST 

INPUT “VARBLE" . . . 


NTHETA 

NRATIO 

NUMRI 

NUMBET 

numlmb 

1 

20 

1 

1 

1 


♦CARDS 4,5,,, ..DATA INPUT FOR NAMELIST "VARBLE" ... 

...DESCRIPTION OF INPUT PARAMETER ARRAYS ... 

NUMBER OF PARAMETER PARAMETER DESCRIPTION 

DATA ARRAY 
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NTHTTA 
NR ATIO 
NUNRI 
NUMBET 
NUMLMB 


TRATIO 

ERATIO 

VRI 

VBETA 

VLMPDA 


PAD PIVOT ANGLE, THETAP/BETA 
PAD SLOPE FACTOR ,EPSILON/HMIN 
PAD INNER RADII, RI 
PAD ANGLE DIMENSION, BETA 
BEARING NUMBER, LAMBDA 


...INPUT FOLLOWS ... 

SVARBLE TRATI0rl.0,VRI=.5,VBETA-45.,VLMB0A=5D., 

ERATIO-. 5,1. ,1.5,2.,Z.5,3. ,3. 5, 4. ,4. 5, 5. ,5.5, 6. ,6. 5, 7. ,7. 5, 8. , 
8.5,9. ,9.5, IG 


SEND 
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APPENDIX D 
FLOW CHARTS 

All FORTRAN symbols used in these flow charts are defined in 
appendix E, 


53 


Executive Program MAIN2 


ENTRY 
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I 












Initialize pad dimensions & 
mesh nodes (variable grid?) 
Reser iteration count KK=0 


ENTRY RREIAX 




Gas film? 


Linearize coefficients Ci._ & D,,„ at 
KKth iteration: ^ 

CCOF = 




Solve finite difference Reynolds Equation at 
KKth iteration by Gauss-Seidel Method with 
relaxation factor for dependent variable: 


Calculate the maximum error between iterations 
ERROR = niaxfU,^^ ~ ^km“^ 


Is 

^x^s (ERROR) 

X^ess than a preset^ 
gst convergence valuej 


Calculate & store pressure array: 






liquid 


RETURN 
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jiquid film 


Calculate edge flow rates 
at R=R,* & R=1 







r A 

RSIMP 

V 7 


Initialize pad dimensions & 
mesh nodes 


R ^ 1 ; O^e ^/3 
i m ’ m r 


Simpson integration along radii; 
parameters Q 


Kir n 

XF = J P. ^ R 
m km 


n=l for pad loading 
2 for moment about 
origin 

3 If 

XF = J (shear stress) dR 
m Ri 


iqui4 film 


4 4 

XFq & XF for edge leakage (flow) 

at leading & trailing edge,resp. 

5 6 • • 

XFjji & XFj^ = radial pressure gradients 

at R=Rj[ & R=l,resp, ,for leakage flows 

at those radii 


7 7 

RETURN 

7 7 
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APPENDIX E 
FORTRAN SYMBOLS 

A(15, 21) coefficient array for U. , . 

AA integral for calculation 

ANG difference 6-6 

P 

AN GB angle /3 in radians 

B(15, 21) coefficient array for U. , . 

1“ > J 

BETA sector dimension (3 

C(15, 21) coefficient array for U. 

CC(5, 3) coefficient array for dP/d6 in fluid shear calculation 

D{15, 21) coefficient array for U. n 

1, J“ J. 

DA A6 of mesh 

DEBUG extensive printout logical switch; TRUE causes detailed output 

of iterative solution routine. FALSE prevents printout. 

DR AR of mesh 

E(15, 21) coefficient array for U. . 

J 

END(4) storage array in numerical integration procedure 

EPS current value of clearance parameter e 

ERATIO(25) input array of e/H^ values 

fit 

ERROR difference between Uv. values at k*'" iteration and (k - 1)_4. 

1] St 

iteration 

EVEN (4) storage array in numerical integration procedure 

EXIT(5) array of integration results from performance calculations 

F(15, 21) array of values for ri^t hand side of nonhomogeneous Reynolds 

equation (liquid film) 

F(K, Ij J) a DEFINE (local) function in radial integration subroutine for 

load (K=l) and moment (K=2) calculation 

FF integral of frictional energy dissipation 
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G(15, 21) 

H(15, 21) 

HALT 

HLMBDA 

HMIN 

HRATIO 

H2(I, J) 

ISTORE 

ITERMX 

JSTORE 

K 

LAMBDA 

LASTA 
LASTR 
LN(I), LS(I) 

LE(J), LW(J) 

M 

NA 

NR 

NRATIO 

NTHETA 

NUMBET 

NUMLMB 


array of factors for nonlinear part of coefficients to 


array of values of dimensionless film thickness h/h^ 

test value for smallest HMIN 

Ah2 

minimum film thickness ratio, h„/h 

maximum- to- minimum film thickness ratio hj^/hg 

a DEFINE fimction in RSIMP for ratio h/hg 

i^^ mesh position where occurs the maximum U. . 
change between successive iterations 

input value for maximum allowed number of iterations 

mesh position where occurs the maximum U.j change 
between successive iterations 

index in DEFINE function F(K, I, J) 

compressibility factor, A; sometimes called the ’’bearing 
number' ' 

trailing edge boundary node 
outer radial boundary node 

finite difference indices for 6P/8R and 8^P/9R^ at R(I) 
for variable mesh coding 

finite difference indices for BP/dO and B^P/dO^ at 0(J); 
for variable mesh coding 

index in DEFINE function X(M, I, J) 

number of angular mesh increments 

number of radial mesh increments 

number of e/Hg values input in ERATIO array 

number of 0p//3 values input in TRATIO array 

number of /3 values input in VBETA array 

number of A values input in VLMBDA array 
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NUMRI 

ODD(4) 

OLDQ 

OMEGA 

P(I, J), PNORM(I, J) 

Q(I, J) 

QSTORE 

R(15) 

RESIDL 

RI 

RR 

STRERR 

TH(21) 

THPR 

TRATIO(IO) 

VARGRD 


VLMBDI(IO) 

VRI(IO) 

WUNIT 


ww 


number of radius ratio values in VRI array 

storage array in numerical integration procedure 

logical switch for use of previous numerical solution 
array of as first guess for new calculation. 
FALSE returns ambient pressure as initial guess. 

iteration relaxation factor, 


normalized lubricant film pressure relative to ambient 
pressure 

the working storage array for the U.j in the compres- 
sible gas program 

value of the maximum change of U.^ for one iteration 
over the entire sector pad domain 

array of values r/r„ at the radial mesh nodes 
I=1(1)LASTR 


input error limit on maximum allowed difference of 
between successive iterations 

current value of radius ratios from VRI array 

integral of first moment for R 

storage value of ERROR during search for maximum 
error within one iteration pass 

array of values for all j-mesh 

nodes J=1(1)LASTA 

input array of 9 ^/^ values 



variable mesh logical switch; TRUE causes finer finite 
difference mesh next to pad boundaries; FALSE utilizes 
a constant mesh increment in either coordinate 

input array of A values 

input array of radius ratio values, RI 

xinit load capability of bearing, i. e. , total load WW per 
unit area 

integral for bearing load capability 
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X(M, I, J) 
XF(21, K) 


a DEFINE (local) function in radial integration subroutine 
RSIMP for friction dissipation. M=1 at J=1,M=2 
for J= 2(1)IASTA-1, and M=3 at J=LASTA 

array of three functions from radial integration subroutine 
with values at each mesh point. K = 1, 2, 3 corre- 
spond to load, center of pressure, and friction calcula- 
tions, respectively. K = 4, 5, 6 correspond to edge 
leakage. 
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APPENDIX P 
SAMPLE PROBLEMS 

Computer output listings are presented for two representative runs, one 
for a compressible lubricant case and one for liquid lubricant. Input data cards 
for both cases are shown with the computer program listings in appendix C. 

In both runs ft=45®, 0Vfl=l. , R,=. 5, and 20 values of e/Hn=. 5(. 5)10. 

pi u 

In addition A=50. for the compressible case. 

The execution of all cases by the compiled programs required less than 
1 minute of computer time on the UNIVAC 1100/42. 


NUMERICAL results - GAS BEARING ANALYSIS 
LEh'Ij RESEARCH CENTER (NASA) 


SYSTEM PARAMETERS 



(PI) 


(P2) 



(P3) 



(P4> 



INNER/OUTER 

RADII 

BEARING 

PAD ANGLE 


BEARING NUMBER 


PIVOT 

angle RATIO 



RI/RO 


BETA (DEGREES) 


lambda 


TKETAP/faETA 



.50000 


45 

.000 


50.000 



1.0000 






NUMERICAL 

RESULTS 






(R1 ) 

(R2) 

(R3) 

(R4) 

(R5) 

(R6) 

(R7) 

(R8) 

(R9) 

(RIO) 

(Rll 

H1/H2 

1/H2 

EPSIL0N/H2 

LOAOtUBAR 

UNIT-LOAD 

FRICTION 

,F F/UBAR 

H2/H(CP) 

R(CP) 

THETA(CP)/ 

X(CP) 










BETA 


1.3536 

1 .0000 

.50000 

.19982-01 

.67847-01 

1.3367 

66.893 

.90113 

.77744 

.63565 

-.21945 

1.7071 

1.0000 

1.0000 

.32872-01 

.11161 

1.2000 

36.504 

.82331 

.77395 

.64226 

-.21461 

2.C607 

1 .0000 

1.5000 

.40893-01 

.13885 

1.0992 

26.879 

.76015 

.77145- 

.64836 

-.21036 

2.‘U‘i2 

1.0000 

2.0000 

.45591-01 

.15479 

1.0214 

22.403 

.70773 

.76960 

.65414 

-.20649 

2.7678 

1 .0000 

2.5000 

.48056-01 

.16317 

.95907 

19.957 

.66346 

.7681» 

.65967 

-.20290 

3.1213 

1.0000 

3. COCO 

.49046-01 

.16652 

.90772 

18.508 

.62555 

.76708 

.66496 

-.19953 

3.<*7<»9 

1.0000 

3.5000 

.49071-01 

.16661 

.86441 

17.616 

.59267 

.76619 

.67000 

-.19637 

3.8284 

1.0000 

4.0000 

,48472-01 

.16458 

.62718 

17.065 . 

.56384 

.76545 

.67480 

-.19339 

4. 1320 

1.0000 

4.5000 

.47480-01 

.16121 

.79466 

16.737 

.53833 

.76484 

.67936 

-.19058 

4. 5355 

1.0000 

5.0000. 

.46247-01 

.15702 

.76589 

16.561 

.51555 

.76431 

.68369 

-.18793 

4.8891 

1.0000 

5.5000 

.44875-01 

.15236 

.74016 

16.494 

.49508 

.76386 

.68779 

-.18545 

5.2426 

1.0000 

6.0000 

.43433-01 

.14747 

.71693 

16.506 

.47655 

.76346 

.69168 

-.18307 

5.5962 

1.0000 

6.5000 

.41968-01 

.14249 

.69560 

16.579 

.45967 

.76311 

.69537 

-.18084 


ORIGINAL PAGE IS 
OF POOR QUALITYi 



5 . 9<»97 

1.0000 . 

7.0000 

. 90508-01 

.13759 

.67695 

16 .699 

.99923 

.76280 

.69888 

-.17873 

6.3033 

1.0000 

7.5000 

. 39075-01 

.13267 

.65862 

16.655 

.93003 

.76252 

.70220 

-.17672 

6. 6569 

1.0000 

8.0000 

. 37680-01 

.12799 

.69210 

17.091 

.91691 

.76227 

.70537 

-.17982 

7.0109 

1.0000 

8.5000 

. 36332-01 

.12336 

.62679 

17.250 

.90976 

.76209 

.70838 

-.17301 

7.3690 

l.OCCO 

9.0000 

. 35036-01 

.11696 

.61239 

17.979 

.39399 

.76183 

.71125 

-.17130 

7. 7175 

1.0000 

9.5000 

. 33793-01 

.11979 

.59695 

17.729 

.38288 

.76169 

.71398 

-.16966 

8.0711 

1.0000 

10.000 

. 32609-01 

.11070 

.58630 

17.983 

.37300 

.76196 

.71659 

-.16810 





NUMERICAL RESULTS - INCOMPRESSIBLE 
LEWIS RESEARCH CENTER ( NASA ) / 

ANALYSIS 

S 
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<P?> 


(P2) 

SYSTEM PARAMETERS 

(P3) 



(?4) ' 

81 
S ^ 

INNER/OUTER 

RADII 

BEARING 

PAD ANGLE 

PIVOT ANGLE RATIO 



& ^ 
a ^ 


RI/RO 


BETA (DEGREES) 


THETAP/BETA 





.50000 

- 

45. 

COO 


1.0000 


• 






'> 

NUMERICAL 

RESULTS 






(Rl) 

(R2» 

(R3) 

(R4» 

(R5) 

(R6) 

( R 7 1 

(R3) 

<R9) 

(RIO) 

(RID 

EPSILON/H2 

H1/H2 

UNIT-LOAD 

FRICTION/ 

RO SIDE 

RI SIDE 

'edge 

INSIDE 

(R(CP)-RI 

)/ THETA (CP)/ 

X)CP) 




LOAD 

leakage 

LEAKAGE 

LEAKAGE 

FLOW 

U.-RI) 

beta 


.5D000 

1.3536 

.23536-02 

40.132 

.43651-01 

.14749-01 

.46783 

.38495 

.54650 

.53403 

-.27671 

l.OODO 

1.7071 

.34440-02 

24.726 

.89588-01 

.29153-01 

.45923 

.39622 

.54107 

.55810 

-.26209 

l.SQOO 

2.0607 

.39366-02 

19.844 

.13708 

.43063-01 

.45196 

.40841 

.53746 

.57776 

-.25029 

2.0000 > 

2.<tl<t2 

.41414-02 

17.534 

.18618 

.56982-01 

.44573 

.42187 

.53474 

.59382 

-.24067 

2.5000 

2.7678 

.42036-02 

16.228 

.2-'734 

.71215-01 

.44030 

.43701 

.53253 

.60688 

-.23285 

3.0000 

3.1213 

.41838-02 

15.448 

.28971 

.85930-01 

.43541 

.45228 

.53057 

.61884 

-.22569 

3.5000 

3. *(749 

.40978-02 

15.056 

.34215 

.10026 

.43113 

.46703 

.52912 

.62988 

-.21913 

M.ooon 

3.8284 

.39756-02 

14.908 

.39492 

.11407 

.42740 

.48190 

.52809 

.63976 

-.21330 

4.5CD0 

4. 1820 

.38353-02 

14.928 

.44785 

.12746 

.42413 

.49672 

.52735 

.64874 

-.20602 

3. GOOD 

4.5355 

.36918-02 

15.052 

.50137 

.14079 

.42122 

.51193 

.52675 

.65683 

-.20327 

5.5303 

4.8691 

.35519-02 

15.248 

.55575 

. 15432 

.41860 

.52777 

.52621 

.66409 

-.19900 

6.0000 

5.2426 

.34119-02 

15.530 

.61004 

.16758 

.41627 

.54337 

.52581 

.67087 

-.19502 

6.5000 

5.5962 

.32794-02 

15.859 

.66515 

.18106 

.41417 

.55955 

.52544 

.67704 

-.19140 

7.0000 

5.9497 

.31492-02 

16.258 

.71996 

.19421 

.41231 

.57536 

.5251<6 

.68287 

-.18798 



7.500D 

6.3033 

.30279-02 

16.690 

.77568 

8.0000 

6.6569 

,29098-02 

17.189 

.83091 

B.SOOO 

7.0X0«» 

.28005-02 

17.703 

.88719 

9.0Q00 

7.3690 

.26949-02 

18.280 

.99280 

9.S0C0 

7.7175 

.25965-02 

18.877 

.99960 

10.000 

8.0711 

.25099-02 

19.502 

1.0571 

aFRCe LIQUIO. 






.20766 

.91061 

.59180 

.52992 

.68819 

-.18986 

.22075 

.90910 

.60779 

.52979 

.69327 

-.18189 

.23921 

.90773 

.62999 

.52955 

.69793 

-.17916 

.29728 

.9P6S1 

.69066 

.52993 

.70291 

-.17653 

.26078 

.90539 

.65760 

.52930 

.70652 

-.17912 

.27960 

.90938 

.67508 

.52916 

.71035 

-.17187 
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Figure 1. - Geometry of sector pad. 





Figure 2. Sector Pad Mesh Definition for Finite 
Difference Scheme. 



LN(I)-=I+1 
LS (I) =1-1 
LE(J)=J+1 
LW(J)=J-1 


(a) constant increment 


LN(I)=I+1 

LS(I)=I-2 

LE(J)=J+1 

LW(J)=J-2 


LECJ) 


(b) variable increment 


1 Derivatives. 




